Boundary Energies and the Geometry of Phase Separation in Double-Exchange 

Magnets 
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We calculate the energy of a boundary between ferro- and antiferromagnetic regions in a phase 
separated double-exchange magnet in two and three dimensions. The orientation dependence of 
this energy can significantly affect the geometry of the phase-separated state in two dimensions, 
changing the droplet shape and possibly stabilizing a striped arrangement within a certain range 
of the model parameters. A similar effect, albeit weaker, is also present in three dimensions. As 
a result, a phase-separated system near the percolation threshold is expected to possess intrinsic 
hysteretic transport properties, relevant in the context of recent experimental findings. 

PACS numbers: 75.47.Gk, 75.30.Kz, 75.10.Lp 



Phase separation is ubiquitous in systems involving 
multiple types of degrees of freedom and opposing in- 
teractions. Among the growing number of compounds 
which are known to exhibit experimental signatures con- 
sistent with this phenomenon one finds the doped man- 
ganese oxides^ [renowned due to their colossal magne- 
toresistance (CMR)], as well as lightly-doped magnetic 
semiconductors^ (e.g., Eu-based). The basic physics of 
both of these systems involves conduction electrons that 
hop between sites in a crystal, whose ionic (localized) core 
spins are coupled to the electronic spins via a strong on- 
site ferromagnetic Hund's exchange interaction. Under 
broad conditions, the energy of the Fermi sea of conduc- 
tion electrons is then minimized when the core spins are 
ordered ferromagnetically, giving rise to double- exchange 
ferromagnetisrrA This ordering is opposed by strong an- 
tiferromagnetic tendencies present in such systems, which 
originate from Heisenberg superexchange between the 
core spins, and which may also be complemented by 
other, band-structure and orbital-related, mechanisms. 
Phase separation resolves this competition by dividing 
the system into regions of different magnetic orders and 
electronic band- widths. We will be concerned with the 
case where the ferromagnetic regions in which the elec- 
trons are free to move (and thereby to lower their kinetic 
energy) co-exist with electron-poor domains where anti- 
ferromagnetism prevails. 

The present article focuses on the properties of the 
boundaries between different regions in the sample. It 
is shown that with increasing strength of antiferromag- 
netism the boundary energy becomes increasingly depen- 
dent on the orientation of the boundary with respect to 
the lattice. Thus, when the Eu chalcogenides (which 
have Neel temperatures of the order of 10 K or below) 
are doped with carriers, the latter are predicted to form 
droplets of circular or spherical shape. In the CMR man- 
ganates, on the other hand, antiferromagnetism is much 
more pronounced and we expect a strong dependence of 
the boundary energy on orientation. This fact, in turn, 
is likely to contribute to the unusual transport proper- 
ties of the CMR manganates near percolation threshold 
including hysteresis, memory and switching effects^. 



Although to the best of our knowledge such considera- 
tions have not been presented before, it appears reason- 
able to expect that at the qualitative level these are not 
specific for the simplified model considered below, nor 
indeed for double-exchange magnetism. Support for this 
view can be drawn from the numerical studies of phase 
separation in the t — J and Hubbard models, showing 
diamond-shaped islands 6 or an anisotropy in the stripe 
direction 7 . Rigorous results on the preferred wall di- 
rection between ferromagnetic domains in the Falicov- 
Kimball model? are likely to be indicative also of the 
behavior of phase separation boundaries in this model. 

We begin with the standard double-exchange Hamilto- 
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Here Cj a (a =T,i) are the electron annihilation opera- 
tors, Si are the operators of the core spins located at the 
sites of a square (or simple cubic) lattice, and the vec- 
tor d a P is composed of Pauli matrices. The ionic spins, 
which are treated classically (assuming S 1), interact 
with each other via a nearest-neighbor superexchange J 
and with the spins of the conduction electrons via a large 
Hund's coupling Jh- We consider the limit 9 Jh — >■ oo 
and omit all peculiarities of lattice and orbital structure 
(which vary between compounds). While the long-range 
Coulomb interaction is absent from the Hamiltonian its 
effects will be discussed qualitatively. Throughout this 
paper we assume zero temperature and take the lattice 
spacing, the bare hopping amplitude t, and h to be 1. 
Since the resulting model is symmetric with respect to 
quarter-filling x = 0.5, we treat the x < 0.5 case only. 

As Jh — ► oo, there is only one spin state available for 
a conduction electron at each lattice site. The carriers 
can therefore be treated as spinless, and their nearest- 
neighbor hopping amplitude depends on the local spin 
configuration via the double-exchange mechanism^. In 
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particular, it equals 1 when the ionic spins at the two 
neighboring sites are parallel to each other, and it van- 
ishes in the case of anti-parallel spins. When a certain 
number of carriers, N e — xN , where N is the number of 
sites in the system, are doped into an insulating (antifer- 
romagnetic) double-exchange magnet, it proves energeti- 
cally favorable for all these carriers to concentrate in one 
part of the sample. The magnetic ordering in this part 
is changed to a ferromagnetic one, while the rest of the 
sample remains antiferromagnetic and free of conduction 
electrons (as long as x and J are not too large, see be- 
low). The optimal value of the carrier concentration in 
the ferromagnetic region, xfm > x, is determined by the 
condition 




where JC and £ are the elliptic integrals^. 

In 3D, the analogue of such an abrupt diagonal bound- 
ary is an abrupt planar boundary perpendicular to a 
main diagonal [direction (1,1,1) or equivalents] of the 
cubic lattice. The calculations of Ref. [Tj can be gener- 
alized to the 3D case, yielding the boundary energy, 
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where in the case of dimensionality D the thermodynamic 
potentials of the antiferromagnetic and ferromagnetic 
phases are given by VLafm{J) = — DJ and Q,fm{J,h) = 
f!^ D (e — n)vD{t)de + DJ. Here the chemical potential 
is denoted by /i — Jfj/2 and ^£)(e) is the density of 
states for the 2D (3D) tight-binding electronic disper- 
sion, e ^ = — cos k x — cos k y (— cos k z ) . 

The values of \x and xfm = J_fi v D{()de depend solely 
on J (see the dotted lines in Fig. [21 below). When the 
value of the nominal density x reaches xfm(J), the entire 
sample turns ferromagnetic. However, an investigation of 
the energy balance and stability of different phases shows 
that thermodynamic equilibrium between Neel and fer- 
romagnetic phases is possible only as long as Xfm is not 
too large, Xfm < x c or equivalently J < J c . While the 
precise value of J c is not known, a variational study^i in 
2D shows that x c < 0.245, corresponding to J c < 0.036. 
In 3D, x c < 0.291 or J c < 0.035, at which point the spin 
stiffness in the ferromagnetic phase turns negative. If the 
value of J exceeds J c the magnetic ordering in either the 
electron-rich or the electron-poor regions of the sample 
will differ from that of a ferromagnet or a Neel antiferro- 
magnet. Rather, one of the many possible antiferromag- 
netic phases characterized by non-vanishing conduction 
electron density may be stabilized (see Ref. [ll| for exam- 
ples). Throughout the present article, we assume that J 
is smaller than J c . 

The presence of phase separation inevitably raises 
questions regarding the structure of the boundary be- 
tween the phases. In an earlier articled, it was suggested 
that the most energetically favorable structure of such a 
boundary corresponds to an abrupt change in the mag- 
netic order. For the case of phase separation into ferro- 
magnetic and Neel phases, which is of interest to us here, 
straight boundaries running in the direction of a lattice 
diagonal in a 2D system were considered [see Fig. ^ ( a )]- 
The energy of such a boundary per unit length is given 
byi£ 



V2Wf D) (J) 



2tt 



m m 

— arccos— 
2tt 2 



\xxfm (3) 



per unit area. In Eq. @, the integration variable k± = 
{hi, k 2 } varies over the first Brillouin zone of a triangular 
lattice of unit spacing. The k\ axis is chosen in one of 
the nearest-neighbor directions. 
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FIG. 1: Diagonal (a) and vertical (b) boundaries in the 2D 



In the case of small J C 1, these somewhat cumber- 
some expressions reduce to 
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In both 2D (Refs. I 111 1 3f) and 3D cases, one can ver- 
ify that tilting the spins adjacent to a diagonal boundary 
(effectively spreading it over two or three lattice links) in- 
creases the boundary energy. This strongly supports the 
notion of an abrupt inter-phase boundary. The fact that 
the boundary is sharp (as opposed to a smooth Bloch-likc 
wall) implies the absence of a continuum limit and hence 
a non-trivial dependence of the boundary energy on ori- 
entation. We will now explore this property variationally 
in more detail. 

Specifically, we consider a boundary running perpen- 
dicularly to the (1,0) and (1,0,0) directions in the 2D 
[FigdfM and 3D cases, respectively. The energy of such 
a "vertical" wall in D — 2, 3 dimensions can be evaluated 
exactly within the same technique (see Appendix). We 
find 

W v = - v D -i(e±)de x i{e,e!_)de-DJ. (6) 
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Assuming that e < 0, the spectral shift function £(e, e±) 
is given by 
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where in the present case Q = 1 an d e± =eT £ i- 
At small J, one obtains 
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The leading-order terms in Eqs. JSJ and © originate 
from the energy of a partition inserted into an ideal gas 
of electrons. In the limit of small J (corresponding to 
small xfm and to a large Fermi wavelength) this quan- 
tity should not depend on the orientation of the partition 
with respect to the lattice, as is manifest from Eqs. |J5J 
and (JHJ. This independence of the boundary energy on 
orientation at small J is expected to hold generally, and 
not only for the two directions considered above. The 
droplets of a ferromagnetic phase within an antiferromag- 
netic sample will then have a circular (spherical) shape, 
obtained by minimizing the boundary length or areaii 
This changes with increasing J as we discuss below. 
The 2D case. The energy of a vertical boundary al- 
ways exceeds that of the diagonal one [see Fig. |21 (a)], 
and with increasing J the ratio w = W v /Wd increases. 
The droplet shape in 2D then approaches that of a dia- 
mond with rounded corners. The curvature radii at the 
corners decrease as w increases. We model this transition 
approximately, allowing only for the vertical and diago- 
nal orientations of the boundary (see Fig. Minimizing 
the boundary energy of the droplet with respect to the 
ratio a of the combined length of the vertical components 
to that of the diagonal ones while keeping the area of the 
droplet constant yields a = (V2 — w)/(V2w — 1). As the 
value of w increases from 1 to y/2 (the latter corresponds 
to J w 0.030 or xfm ~ 0.221), the droplet shape changes 
from that of a regular octagon^ to that of a diamond. 

This conclusion is expected to affect the properties of 
the system in a case where the ferromagnetic phase oc- 
cupies a considerable fraction of the sample. The long- 
range Coulomb interaction dictates that instead of a 
single charged macroscopic droplet of a ferromagnet, a 
Wigner crystal-like array of smaller droplets is formed^. 
Provided that these droplets are large enough^ and well 
separated in space, the shape of each droplet is governed 
by the results reported above. The small-J regime of 
w — 1 <C 1, and hence circular droplets, is predicted for 
Eu-based semiconductors^. As more carriers are doped 
into such a system, the geometry of the state eventually 
changes from that of disconnected droplets of a (metallic) 
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FIG. 2: The dependence of the diagonal, vertical, and [in (b)] 
{110} boundary energies (solid, dashed, and dashed-dotted 
lines, respectively) on the value of the superexchange coupling 
J in the 2D (a) and 3D (b) cases. Dotted lines (right scale) 
show the carrier density Xfm(J) m the ferromagnetic regions. 
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FIG. 3: Evolution of the droplet shape in the 2D case. 



ferromagnet embedded in an antiferromagnetic insulator 
to the inverse situation of a connected metallic phase 1 . 

The regime of larger J > 0.015 is expected to be re- 
alized in the CMR manganateait The corresponding w 
exceeds 1.25, leading to the following consequences: 

(i) Hysteretic effects in transport. When the ferromag- 
netic droplets are not connected an applied voltage gives 
rise to internal electrostatic fields. Connected conduction 
paths could then form in order to lower the associated 
field energy. When subsequently the voltage is applied 
along another direction, these paths should restructure 
accordingly. Due to the dependence of the boundary en- 
ergy on orientation, this cannot occur smoothly, but nec- 
essarily involves discrete transitions (over energy barri- 
ers) and hysteresis. Experimentally, switching and mem- 
ory effects in the resistivity of CMR manganates near the 
percolation threshold have been observed both in thin 
films^ and in 3D single crystals^. We emphasize that 
such properties arise here as intrinsic to the system, as 
opposed to more conventional effects of impurity pinning 
of the boundary. 

(ii) For a fixed area fraction of the ferromagnet, x/xfm, 
the ratio A = (Cl s — ildr)V^/ {e^Wd) 1 ' 2 , where fldr and fl s 
are the energies of a droplet phase and of a phase com- 
posed of diagonal ferro- and antiferromagnetic stripes, 
depends solely on the droplet geometry! 2 , (here, e is 
the dielectric constant). Variational studies!*^ suggest 
that for w = 1 (circular droplets) and x/xfm ~ 1/2 
the quantity A is only slightly larger than 0. While 
leaving this issue for future study, we note that since 
^drVe/(e 2 Wd) 1 ^ 2 is expected to be larger for diamond- 
shaped droplets (due to a higher value of boundary en- 
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ergy per particle) , it is possible that for large w a diagonal 
stripe phase proves energetically favorable (A < 0) in a 
region near x/xfm — 1/2. In a broader physical context, 
droplet shape effects may complement other mechanisms 
suggested^ to stabilize a stripe phase in the CMR man- 
ganates. 

The 3D case. As seen from Fig. [5] (b). the energies 
of the vertical and diagonal boundaries follow each other 
closely (with W d > W v for J < 0.00019) throughout the 
relevant range of the values of J. Thus, using the same 
arguments as above (cf. Fig. |3J) we conclude that the 
optimal droplet shape remains close to that of an octa- 
hedron (formed by diagonal boundaries) with its vertices 
"shaved off" by vertical and horizontal planes. Energeti- 
cally, this geometry is nearly indistinguishable from that 
of a sphere. In such a case, one can expect based on 
the variational results for spherical droplets^*^ that an 
alternating arrangement of ferro- and antifcrromagnetic 
planar slabs is unstable for all values ofx/x fm ■ However, 
the boundary energy even in 3D is not at all orientation- 
independent. This is illustrated by the dashed-dotted 
line in Fig. [2] (b), which shows the energy Who of a 
boundary perpendicular to the (1,1,0) directional, 
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Here, the function £ is given by Eq. @ with Q = 
2 cos(fcj,/2) and e± = e ± cos k z . At small J, we find 
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As illustrated by Fig. |21 (b), Wxio exceeds both W4 and 
W v by some 20% in an extended range of values of J. 
Thus, our conclusion on hysteresis due to abrupt switch- 
ing of the conduction paths persists in 3D (at least for 
some directions of the path). This has an enhanced sig- 
nificance because in 3D, point defects are less likely to 
interfere with the restructuring of conduction paths. 

We note that in a phase-separated system quenched 
(chemical) disorder that leads to variations in the on-site 
energies also tends to sub-divide large ferromagnetic ar- 
eas into droplets. Since the boundary energy originates 
from an integral of the carrier phase shift over the entire 
Fermi sea, it is expected to be only weakly sensitive to the 
detailed structure of individual wavefunctions (as long as 
the disorder is not too strong). Hence, we expect the ori- 
cntational dependence of the boundary energy to persist 
in the presence of such inhomogeneities (and of many- 
body effects due to the short-range part of the Coulomb 
interaction^) . This is illustrated by numerical-simulation 
results for a 2D phase-separated double-exchange mag- 
net with disordei— , showing that the prevailing direction 
of the boundaries is clearly along the diagonals. 

We note that in order to make quantitative compari- 
son with the experimental results for the manganates one 
may have to complement Eq. with further terms re- 



flecting orbital and lattice effects inherent to these com- 
pounds. While these may affect the optimal shape or 
the preferred direction of the boundaries in a phase sep- 
arated system, our overall conclusion on the presence of 
an orientational dependence of the boundary energy is ex- 
pected to persist. As demonstrated above for the simple 
model this general conclusion follows from the pre- 
sumed short-wavelength (abrupt) nature of the bound- 
ary, combined with sufficiently large band-filling. The 
latter ensures that the short-wavelength (on the scale of 
the lattice constant) properties of the electron gas are 
sensitive to the orientation with respect to the lattice. 
Including phonon or orbital degrees of freedom may mod- 
ify this short wavelength behavior, but is rather unlikely 
to eliminate it altogether. 

As we have previously mentioned, experimental ob- 
servations of non-ohmic behavior and memory effects in 
the resistivity of CMR manganates near the percolation 
threshold exist both for 3D single crystals^ and for thin 
films 5 ". While qualitatively this is in line with our expec- 
tations, it would be desirable to have direct imaging of 
the spatial structure of the underlying phase-separated 
state (as provided by magnetic force microscopy). Such 
measurements could also be used to search for any puta- 
tive stripe phases, especially in 2D. We are not aware of 
similar transport studies of the Eu-based magnetic semi- 
conductors, and indeed we do not expect intrinsic non- 
linear effects to be prominent in this case. However, this 
point should also be confirmed by experiments. 

We acknowledge helpful discussions with O. Agam, G. 
Jung, and V. Markovich, as well as the support of the 
ISF (under the Centers of Excellence Program) , the BSF 
(grant No. 2004162), and the Israeli Absorption Ministry. 



APPENDIX 

In the following we briefly outline the calculation of 
the "vertical" boundary energy, Eqs. (JSH7|) using the 
theory of local perturbations due to I. M. Lifshits and 
M. G. Krein22. We first evaluate the double-exchange 
(electronic kinetic energy) contribution to the boundary 
energy. The (auxiliary) unperturbed system consists of 
two disconnected components: 

(i) A square (in 3D, cubic) tight-binding network of N 
sites, which we divide into two sub-lattices, labeled b and 
c (see Fig. |3J, thereby doubling the period in the direc- 
tion perpendicular to the x axis (in 3D, this corresponds 
to a checkerboard arrangement in the yz plane). Per- 
forming a partial Fourier transform in the y (in 3D, y 
and z) direction, renders the Hamiltonian as 



-\ h \S x )hS x + !) - \ c \S x ^S x + V 



* = E l _ 

x,kj_ 

+e ± (k ± )ct (x)bz (x)+H.c 



(A.l) 
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where the operators b and c annihilate the two species of 
spinless fermions on the two sub-lattices, and e±(k±) = 
— cos k y {— cosfcz) is the usual tight-binding dispersion in 
D — 1 dimensions. 

(ii) A chain of \/N sites (for the 3D case, a plane of 
N 2/>3 sites), disconnected from each other (zero hopping 
amplitudes, hence H2 = 0). We again introduce two 
sub-lattices in a checkerboard fashion, and perform the 
Fourier transform, labeling the two resulting localized 
fermion operators / and g with the D — 1-dimensional 
momentum index, fcj_. 
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FIG. 4: The effect of the perturbation, Eq. E"2t . in the 2D 
case. Carrier hopping along the lattice links marked "x" is 
eliminated, and the new bonds (dashed lines) are established. 

Next, consider a perturbation of the form 



V 



The first two terms correspond to disconnecting the x = 
chain (plane) from the rest of the system. An extra site 
is then connected to every other site at x = ±1 (last 
two terms in V), creating a vertical boundary perpen- 
dicular to the x axis (see Fig. 0J. Note that hopping 
is allowed both in the right, and in the left half-plane 



(half-space). Therefore, the change in energy of the sys- 
tem due to the perturbation (when supplemented with 
the appropriate superexchange contribution) is equal to 
twice the energy of a vertical boundary between an anti- 
ferromagnet (where the hopping amplitude vanishes) and 
a ferromagnet, as pictured in Fig. ^ (b)- 

Since the transverse momentum k± is preserved by V, 
the variables separate. For a given k±, the perturbation 
is localized at x = 0, ±1 and takes the form 



i=l 



id] a 



(A.3) 



where 2A,- t = { — -\/3, —1,1, V3, — \/2, V2} and a, are suit- 
able combinations of the original fermion operators that 



satisfy a^at + a\ai 



The spectral shift function, Eq. 



0, is then evaluated as 22 
1 



with 
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f(fci,e) = -^Arg Det {l - M, r \,) , 
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Here, the Dirac bra and ket states are defined by |/g ) = 
ft |0) (with |0), the vacuum state), etc., and the opera- 

rex 

tors pi are the Fourier components (with respect to k x ) 
of (c£ (x) ± 6g (x))/\/2- As before e± = e =F ej_(fcj_). 

The choice of the Arg branch in Eq. I|A.4|) is sub- 
ject to the requirement that t;(kj_,e) — > for a van- 
ishing perturbation strength, e.g., if the Hamiltonian is 
given by Hi + H2 + ctV with a — > 0. In the 3D case, 
we note the appearance of bound (Tamm) surface elec- 
tron states at the vertical border. These have the en- 
ergy E — ±cosh(fca;) + ej_ = with |ej_(A-j_)| > 1. For 
example the wave function of an ej_ > 1 state in the 
x < half of the system differs from zero only on the / 
sites (x = 0), on the c sites with even x, and on the b 
sites with odd x (see Fig. Q), where it is proportional to 
exp[k x x+ik y y+ik z z]. The net spectral shift function, Eq. 
(|A.4p . contains contributions from these boundary states, 
from the extended states in the two half systems and also 
from the states in the disconnected x — chain/plane 
(see Fig. 0J. The latter needs to be subtracted when 
calculating the desired boundary energy. The resulting 
£(e, e±) is given by Eq. 0. 

By counting the number of antiferromagnetic links and 
the associated energy change, one arrives at the final ex- 
pression, Eq. ©. The energy of a {110} boundary, Eq. 
@, can be evaluated in a similar way. 
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